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Abstract: Reducing the acquisition time for two-dimensional nuclear magnetic resonance 
(2D NMR) spectra is important. One way to achieve this goal is reducing the acquired data. 
In this paper, within the framework of compressed sensing, we proposed to undersample 
the data in the indirect dimension for a type of self-sparse 2D NMR spectra, that is, only a 
few meaningful spectral peaks occupy partial locations, while the rest of locations have 
very small or even no peaks. The spectrum is reconstructed by enforcing its sparsity in an 
identity matrix domain with £ p (p = 0.5) norm optimization algorithm. Both theoretical 
analysis and simulation results show that the proposed method can reduce the 
reconstruction errors compared with the wavelet-based l\ norm optimization. 

Keywords: NMR; spectral reconstruction; sparsity; undersampling; compressed sensing 
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1. Introduction 

Nuclear magnetic resonance (NMR) spectroscopy is widely utilized to analyze the structures 
of chemicals and proteins. Multidimensional NMR spectra can provide more information than 
one-dimensional (ID) NMR spectra. The acquisition time for a conventional two-dimensional (2D) 
NMR spectrum is mostly determined by the number of t\ increments in the indirect dimension. One 
possible way is to reduce the acquisition time is to reduce the number of t\ increments. However, this 
will result in aliasing of the spectrum in the indirect dimension [1,2], because the sampling rate is 
lower than the requirement of the Nyquist sampling rule. 

Researchers have been seeking ways to suppress the aliasing from the aspects of sampling and 
reconstruction. Radial sampling presents relatively small leakage artifacts [3] and Poisson disk sampling 
is observed to provide a large low-artifact area in the signal vicinity [4]. The maximum sampling time 
for multi-dimensional NMR experiments was analyzed by Vosegaard and co-workers [5]. Besides the 
sampling patterns, some reconstruction algorithms have been employed to improve spectral quality, 
including maximum entropy [6,7], iterative CLEAN algorithm [8] and Bayesian reconstruction [9]. 
The sparse sampling was incorporated with intermolecular multiple-quantum coherences for 
high-resolution 2D NMR spectra in inhomogeneous fields [10]. 

Recently compressed sensing (CS) theory [1 1,12], for reconstructing signals from fewer numbers of 
measurements than the number that the Nyquist sampling rule requires has attracted lots of attention in 
medical imaging [13], single pixel imaging [14], and computer vision [15], etc. Under the assumption 
that the acquired data is sparse or compressible in a certain sparsifying transform domain, CS can 
successfully recover the original signal from a small number of linear projections with little or no loss 
of information. The choice of sparsifying transform is important in the CS. The sparsfying transform 
should be maximally incoherent with the measurement operator. Intuitively, the target signal should be 
sparsely represented in the transform domain, e.g., wavelet transform domain, and this spare 
representation should be spread out in the encoding scheme. Iddo introduced CS to reconstruct 
a 2D NMR spectrum from partial random measurements of its time domain signal under the 
assumption that the spectrum is sparse in the wavelet domain [16]. 

In this paper, we focus on the reconstruction of self-sparse NMR spectra, that is, a few meaningful 
spectral peaks occupy partial locations while the rest locations have very small or even no meaningful 
peaks. NMR spectra includes regions where no signals arise because of the discrete nature of chemical 
groups [17]. The reason we pay attention to self-sparse NMR spectra is that many NMR spectra of 
chemical substances fall in this type [3,10,16,17]. Based on the concept of sparsity and coherence in 
CS, we demonstrate that a wavelet transform is not necessary to sparsify the self-sparse NMR spectra 
or even worsens the reconstruction. We propose to reconstruct the NMR spectrum by enforcing its 
sparsity in an identity matrix domain with a € p (p = 0.5) norm optimization algorithm. Simulation 
results show that the proposed method can reduce the reconstruction errors compared with the 
wavelet-based l\ norm optimization. 

Recently, Kazimierczuk and Orekhov [18] and Holland et al [19] independently proposed to use 
CS in proton NMR and showed promising results in reducing acquired data. A combination of 
spatially encoding the indirect domain information and CS was proposed by Shrot and Frydman [20]. 
The spectra were considered to be sparse themselves [18-20], differing from the sparse representation 
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using wavelets [16]. However, no comparison on the reconstructed spectra with and without wavelet 
transform was given and no theoretical analysis was presented. In this paper, we will analyze the 
performance of wavelet transform in the CS-NMR basing on the sparsity and coherence properties and 
simulated results. 

The remainder of this paper is organized as follows. In Section 2, the reason to undersample the 
indirect dimension is given by calculating the acquisition time for a 2D NMR spectrum. In Section 3, 
the two key factors of CS, sparsity and coherence, are briefly summarized and their values are 
estimated for 2D spectra, followed by the proposed reconstruction method. In Section 4, reconstruction 
of self-sparse NMR spectra is simulated to show the shortcomings of the wavelet and the advantage of 
the identity matrix. The improvement of utilizing the € p norm is also demonstrated. Finally, discussions 
and conclusions are given in Section 5. 

2. Undersampling in the Indirect Dimension of 2D NMR 

In NMR spectroscopy, a typical sampled noiseless time domain signal can be described as a sum of 
exponentially decaying sinusoids: 



where J is the number of sinusoids, Aj 9 0y, r 7 and coj are the amplitude, phase in radians, decay time and 
frequency, respectively, of the yth sinusoid [21]. At is the sampling interval and k (k = 0, 1, . . ., K) is an 
integer to denote the Ml sample point. Such a signal will give rise to a spectrum that is the sum of 
Lorentizian peaks centered at different frequencies coj [21], where j corresponds to y'th type of nuclear 
spins. A conventional ID single pulse NMR experiment enforces an excitation pulse on a sample 
followed immediately by data acquisition. The signal eventually decays due to relaxation [22], thus it 
is called free induction decay (FID). Fourier transform (FT) is applied on the FID to obtain a frequency 
domain spectrum. Figure 1 shows the simulated FID signal and the corresponding ID NMR spectrum 
obtained from FT. 

Figure 1. Simulated FID data in time domain (a) and its corresponding ID NMR 
spectrum (b). Note: the FID is simulated according to Equation (1) with J = 2, A\ = 0.5, 
A 2 = 1, At = 0.01 s, n = r 2 = 800, 0 ± = 0 2 = 0, and coi = 70 Hz, co 2 = 20 Hz. 
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The typical experimental time for a ID NMR spectrum usually takes several seconds, thus it is not 
time consuming. However, for a 2D NMR spectrum, the time domain signal is generated based on two 
time variables t\ and h. As shown in Figure 2, one scan of 2D NMR spectrum contains three steps: 
first, the sample is excited by one or more pulses in the preparation period. These pulses result 
in the evolution of magnetization with time t\\ then, the sample is further excited in the mixing 
period; finally, an FID signal is recorded as a function of ti. Usually, t\ is set as t\ = At\ 9 2At\ 9 ... 9 n\At\ 9 
N\At\ (The increment At\ is usually at the order of milliseconds). The number of t\ increments (N\) is 
determined by: 



4/T 



(2) 



where SW 1 = — is the desired spectral width and Af ± = — - — is the corresponding spectral resolution. 

The typical N\ is from 50 to 500 [22]. Given a fixed t\ = n\At\ 9 one scan is performed and the FID 
signal is recorded and stored along the direct dimension. After the scan, the nuclear spins are allowed 
to return to their equilibrium states before the next scan for t\ = (n\ + \)At\ [22]. 

Figure 2. General scheme for 2D NMR spectra. 
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Finally, 2D FT is performed on the 2D FID data. If the time for performing all the pulses in one 
scan is t p , the total scanning time for a 2D NMR spectrum will be: 

N, f (i , Ar \ a. A 



T Nl = £ ( d x + n x At, +t m +t 2 +t p ) = N x 



^ L - L + t m +t 2 +t p 



(3) 



In order to obtain a good resolution in the indirection dimension, N\ is usually several tens or 
hundreds or even more. This will cause the total scanning time for a 2D NMR spectrum to be tens of 
minutes or even several hours [22-26]. 

In this paper, we aim to reduce the scan number for the t\ dimension. Rather than using the uniform 
increment in the indirect dimension (t\ = At\, 2At\ 9 n\At\ 9 N\At\) 9 we randomly choose unduplicated 
Q numbers from n q E {1, 2, N\} 9 and let t\ = n q At\. Let: 

Q 

(4) 

be the sampling rate in this paper, the total time to scan a 2D NMR spectrum is approximately: 

T q=§- T n 1 =PT Ni (5) 

The approximation is made by ignoring the total evolution time £n q G {1,2,...,^}^ = 1,2, ...,q n q^i since 
this value is only in the order of seconds. Compared to the time to acquire a 2D spectrum with fully 
sampled FIDs in the indirection dimension, undersampling the FIDs in the indirect dimension can 
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greatly reduce the acquisition time for a 2D NMR spectrum if p is small enough. Figure 3 shows an 
example where we randomly undersample the indirect dimension with sampling rate p = 5/1 1 = 0.45. It 
means we save nearly half of the acquisition time of the conventional scheme. 

Figure 3. An example of random undersampling in the indirect dimension. The symbol <= 
denotes the acquired FIDs. 




However, this undersampling will result in aliasing artifacts [1,6]. It would be of great value if we 
can minimize these artifacts and reconstruct the full 2D NMR spectrum from the limited data. Here we 
explore the undersampling and reconstruction methods under the framework of CS. 

3. Reconstruction of 2D Self-Sparse NMR Spectra with Compressed Sensing 

3.1. Basic Concepts in Compressed Sensing 

The CS proposed by Candes et al. [11] and Donoho [12] is a new theory to do undersampling and 
reconstruct the signal of interest from limited physically acquired data. They build a theoretical 
foundation that one can exactly or approximately recover signals from highly incomplete 
measurements. The two basic tenets to guarantee the performance of CS are sparsity and incoherence. 

(a) Sparsity. For the signal xER N and a basis dictionary VGR SxJV (e.g., identity matrix, FT, 
discrete cosine transform or wavelet transform matrix), the sparsity is often interpreted as: 

HMIo=Mo«* (6) 

where ||a|| 0 denotes the 4 norm that counts the nonzero entries in a, and S is the number of nonzero 
entries. If x is sparse without transformation (namely sparse in identity matrix I G R Nx N ), it is called 
self-sparse since other complicated sparsifying transform, e.g., wavelet transform, is not required. 

Candes et al. [11] and Donoho [12] proved that it is possible to recover the original signal x from 
O(iVlogS) measurements. This means the required number of measurements is proportional to the 
number of nonzero entries in the basis *F. The smaller the S is, the less the number of measurements 
is required. 

(b) Incoherence. When a signal x is sampled by a sensing matrix Omxtv, the measurements y G R M 
of x is: 

f = Ox (7) 

The coherence is defined as [27,28]: 
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^(0,Y) = max|(^,i«r.)| (8) 

where 0 k is the Mi rows of <D and Wj is the yth column of *F. The coherence measures the largest 
correlation between any row of <D and column of *F. The less the coherence between <D and *F is, the 
smaller the // is. The value range of ju is [1, VN]. The minimal coherence ju = 1 occurs when <D and *F 
is a time-frequency pair [29]. CS requires the coherence to be as small as possible, which means each 
measurement vector 0 k must be 'spread out' in the *P domain [28]. 
If the signal x satisfies [30]: 



II Ho 2 



V 



//(0,<F) 



(9) 



it can be perfectly recovered by solving: 

d = mm||a|| o , s.t. y=OTa (10) 

where ||a|| 0 denotes the £o norm that counts the nonzero entries in a. 
The recovered signal is: 

x = x Va (11) 

Equation (9) implies that if the coherence between <D and *P is small, more non-zeros can be 
allowed in the sparse representation a. CS suggests <D to be random enough to guarantee its 
incoherence with any *F. This is also observed that random sampling in time domain can improve the 
quality of reconstructed spectra [31]. 

However, £o norm is known to be intractable and sensitive to noise [11,12], and l\ norm convex 
optimization is commonly used in CS to recover x by solving: 

d = min||a|| l , s.t. y=OTa (12) 

The accuracy of CS reconstruction using Equation (12) can be guaranteed if <D*F satisfies the 
appropriate restricted isometry properties [32]. A restricted isometry constant a s [32] defined as the 
smallest number such that: 

(l-o- s )|a£<|OTa|<(l + o- 5 )||a|; (13) 

holds for all vectors that have at most S nonzero entries. If a 2S < V2 — 1, the solution to the t\ norm 
problem is that of the £o problem [32]. 

The number of measurements M should satisfy: 

M>C'jU 2 (<S>^)S\ogN (14) 

so that the signal x can be exactly recovered from measurements y in overwhelming majority of 
cases [28]. Equation (14) implies that the number of measurements is proportional to the number of 
nonzero entries S in a and the square of coherence ju. If both S and ju are small, the required number of 
measurements M could be small. This means that one can perform fewer measurements to save 
acquisition time while reconstruct original signal x very well. 
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Iddo [16] applied CS to remove the aliasing artifacts from incompletely acquired FID data by 
enforcing the sparsity of 2D NMR spectra in wavelet domain according to: 

& = min|a| 1? s.t. ||y-0F r T r a|| < a 

where y is the measurements in time domain, 0 is a random sampling operator defining the FIDs 
acquired in the indirect dimension, ¥ T denotes the inverse 2D FT, and T r is the inverse 2D wavelet 
transform. According to Equation (11), the recovered spectrum is x = H^a. 

In this paper, we focus on the reconstruction of self-sparse NMR spectra in which significant peaks 
take up partial locations of the full NMR spectra while the remaining locations have very small or even 
no peaks. Ideally, if the number of sinusoids Jin Equation (1) is very small, and the meaningful peaks 
are narrow enough relative to the whole 2D frequency coverage, the spectra can be considered to be 
sparse since the number of non-zeros for the spectra is much smaller than the number of spectrum 
points in the 2D NMR spectra. 

The sparsifying transform and the coherence between *P and <D = 0F r play important roles in the 
CS, as we have discussed. In the following sections, we will demonstrate that wavelet is not necessary 
to sparsify or even worsens the self-sparse NMR spectra based on the concept of sparsity and 
coherence. We will then reconstruct the NMR spectrum by enforcing its sparsity in an identity matrix 
domain with t p (p = 0.5) norm optimization algorithm. 

To represent the NMR spectra in conventional way [4-7,17], the X and Y coordinate axes are shown 
with unit of parts per million (ppm) [21] defined as: 

where 5 is the chemical shift of a peak with frequency a), coref is the frequency of a reference peak and 
coo is the spectrometer carrier frequency. 

3.2. Sparsity of Self-Sparse NMR Spectra 

Figure 4(a) shows a 2D correlation spectroscopy (COSY) spectrum where most of the peaks 

fill partial and very limited regions of the full spectrum. This leads to the sparsity of spectrum because 
the number of non zeros in the 2D spectrum is much smaller than the number of spectrum points. This 
phenomenon is also observed by Yoh Matsuki et al. [17]. 

To test the sparsity of NMR spectra, we can measure the decay of coefficients in a sparsifying 
transform domain and evaluate the approximation error by retaining the £-term largest coefficients, 
because the reconstruction error is proportional to the power law decay k~ r , where r is a constant 
implying the sparsity of signal [29]. Rapid decay of coefficients implies that one can use less non-zero 
coefficients to approximate a NMR spectrum. If we directly measure the decay of signal without 
complicated sparsifying transform, e.g., wavelets, it means measure the self-sparsity of signal. 
Mathematical saying is measuring its sparsity in the identity matrix. 



Sensors 2011, 11 



8895 



Figure 4. Sparsity of a COSY spectrum and its wavelet (symmlet wavelet 

with four decomposition levels and eight vanishing moments) representation, (a) The fully 
sampled NMR spectrum; (b) decay of real part of spectrum and its wavelet coefficients; 
(c,e) reconstructed spectra from 3% and 1% largest coefficients in wavelet domain; 
(d,f) reconstructed spectra from 3% and 1% largest coefficients in identity matrix domain. 
Note: the wavelet fails to represent peaks marked with arrows in (e) and these peaks are 
successfully represented in (f). 
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As shown in Figure 4(b), both the spectra and its wavelet coefficients can achieve rapid decay. By 
retaining 3% largest magnitude coefficients, the spectra can be reconstructed well in Figure 4(c,d). 
However, the spectrum is sparser than its representation in the wavelet domain. This is demonstrated 
by the faster decay of spectrum than that of its wavelet coefficients in Figure 4(b). By retaining the 1% 
largest magnitude coefficients, the wavelet fails to represent some peaks while the spectrum itself can 
represent these peaks, as marked by the arrows in Figure 4(e,f). 

1 13 

For a 2D H- C COSY spectrum, the spectrum decays faster than its wavelet coefficients 
(Figure 5(b)). This implies that the identity matrix can provide a sparser representation of spectra than 
a wavelet does. Peaks are lost or distorted by using the wavelet transform to represent the spectrum 
(Figure 5(e)), but the spectrum is represented very well with the identity matrix (Figure 5(f)). This 
phenomenon is consistent with the observation on the 2D l H- l H COSY spectrum discussed above. 

Figure 5. Sparsity of a ^-^C COSY spectrum and its wavelet (symmlet wavelet 
with four decomposition levels and eight vanishing moments) representation, (a) The fully 
sampled NMR spectrum; (b) decay of real part of spectrum and its wavelet coefficients; 
(c,e) reconstructed spectra from 1% and 0.1% largest coefficients in wavelet domain; 
(d,f) reconstructed spectra from 1% and 0.1% largest coefficients in identity matrix 
domain. Note: the wavelet fails to represent peaks marked with arrows in (e) and these 
peaks are successfully represented in (f). 
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Figure 5. Cont. 
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As a result, this spectrum is self-sparse, which means spectrum is sparse in the identity matrix. Thus, 
according to Equations (9) and (14), it is better to use an identity matrix than to use a wavelet to 
reconstruct the self-sparse spectra from undersampled FIDs since the wavelet cannot provide a sparser 
representation of the spectrum. In fact, Stern et al. [33] proposed to do iterative soft thresholding on 
the spectrum directly, not on wavelet coefficients, to recover one dimensional NMR spectra from the 
truncated FIDs. Although the sparsity of NMR spectra is not explicitly expressed in that work [33], the 
recovered spectrum is obtained from minimizing l\ norm of spectrum, which implies enforcing the 
sparsity of the spectrum. The problem of their method is that truncation violates the random sampling 
scheme in CS and results in strong Gibbs ringing which is hard to suppress [29]. What is more, 
truncating the ID FID is not necessary to save the time to scan a spectrum since scanning a ID NMR 
spectrum is fast and only takes on the order of seconds. 

3.3. Coherence Property of Wavelet-Based and Identity Matrix-Based CS-NMR Spectra 

Besides the sparsity of signal, another key factor for CS is the coherence between <D and *F 
According to Equations (9) and (14), fewer measurements are required for signal sampling system <D if 
it is less coherent with *F and the signal has same sparsity for different *P. 

Pioneering work on CS has pointed out that the coherence of a time-frequency pair is 
ju(<S>, I) = //(0F r , I) = 1 [28]. Thus, we only need to compute the coherence between undersampled 
Fourier operator <t> and wavelet basis *F r . 

The undersampling of 0 in the indirect dimension is carried out by choosing some of the FID points 
in this dimension. To make this undersampling intuitive, a binary mask which has the same size of 2D 
FID is shown as the undersampling pattern in Figure 6(a). If the value of mask at location (ij) is equal 
to 1 shown as a white pixel, the FID at location (ij) is acquired. 

To avoid the influence of randomness on the coherence calculation, 0 is randomly generated 10 times 
and the coherence is averaged for each sampling rate. Figure 6(b) shows that the coherence between 
wavelet and undersampled Fourier operator <D is larger than the coherence between identity matrix and 
<D. So, from the aspect of coherence, it is also better to choose the identity matrix for self-sparse 
NMR spectra. 
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Figure 6. Coherence of wavelet and FT. (a) One sampling pattern in the indirect dimension 
with sampling rate p = 0.30 (fully sampled points in the indirect dimension 
is N\ = 64); (b) coherences for different sampling rates. The symmlet wavelet 
with four decomposition levels and eight vanishing moments is chosen as a typical wavelet 
for test, which is also the typical wavelet in [16]. Error bar stands for the standard deviation 
when repeating 10 times at each sampling rate. 
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3.4. Reconstruction of Self-Sparse NMR Spectra with £ p Norm Minimization 

In this paper, we propose to reconstruct the self-sparse 2D NMR spectra with identity matrix I 
as follows: 



x = min ||x|| , s.t. y= Q>x 



(17) 



where <D = 0F r . 

To further improve the reconstruction, a l p (0 < p < 1) norm is incorporated which has been 
demonstrated to give better reconstruction of MR images with fewer measurements than l\ norm 
does [34-37]: 



x = min x ,s.t. y=Ox 



(18) 



where ||x||^ = 2n=i \ x n\ p an d x n is the nth entry of vector x. For the function X*) = \xf, with p — > 0, 
fx) gets closer to the £o norm of x, as shown in Figure 7. 

Figure 7. The value of fx) = \xf versus the value of p. 
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Theoretically, the required number of measurements [38] by enforcing the sparsity with a l p 
(0 <p < 1) norm is: 

M > Q (p)K + pC 2 (p)K log(7V IK) {l9) 

where C\ and C2 are determined explicitly and bounded in p and the recommend p is 0.5 [34]. 

In this paper, the t p norm minimization is solved via the ^-shrinkage operator [39] with continuation 
algorithm [40] because of its fast computation. This algorithm is abbreviated as PSOCA and 
summarized in Algorithm 1 . 

Algorithm 1. Self-sparse NMR spectra reconstruction with undersampled data using PSOCA. 

Initialization: 

Input the sampled FID data y, set the regularization 
parameter X =10 and tolerance of inner loop r/ = 5 x 10 . 
Initialize x = F0 r y, Xi ast = x, /? = 2 6 , and a = 0. 
Main: 

While p<2 16 
Inner loop: 

1 . Given x, 

For j = 1 to J, solve Equation (20), 
the solution is a; 

2. Given a, 

solve Equation (22), the solution is x ; 

3. If ||Ax|| = 1 1 xiast - x|| > rj 9 xi as t <- x, go to step 1; 

Otherwise, go to step 4; 
Outer loop: 

4. x <- x, /? <— ip, go to step 1 . 
End While 

Output: x 



For a given continuation parameter /?, PSOCA is implemented to solve two sub-problems: 

(1) ^-shrinkage operator 

a j = S p e (x . ) = max {x . - e |x ; \ P ' X , o} ^ (20) 

I A 

1 

where ^ = (3v~ 2 and y5 is a parameter to be updated in the continuation scheme, x 7 and a 7 are the y'th 
entry of column vectors x and a, respectively. 

(2) solve the linear equation: 




which can be simplified to: 

(01 + XP) F r x = /3F T a + Wy 



Sensors 2011, 11 



8900 



where the term P = 0 T 0 is a diagonal matrix consisting of ones and zeros. The diagonal entries of P 



correspond to the location of FID data and the entry value is 1 if a corresponding FID data point is 
sampled, otherwise the entry value is 0. Equation (22) can be solved fast since only a discrete Fourier 
transform and entry-wise division are required. 

4. Simulation Results and Analysis 

In this section, we will show the advantages of the proposed method in two aspects: (1) identity 
matrix as the sparsifying transform is compared with wavelet transform; (2) t p norm minimization is 
compared with t\ norm minimization. The recommended value of p is 0.5 for stability from empirical 
experiments [34]. The notation 4.5 is short for l p with p = 0.5. The typical l\ norm minimization 
algorithms compared in this paper include iterative soft thresholding (1ST) algorithm [16,41-43], 
alternating and continuation algorithm (AC A) [40]. The ACA is just p = 1 in PSOCA. 

Because regions of small spectrum values usually contain no peaks for practical analysis, we set 
magnitude smaller than a constant Tto be zero according to: 



where x denotes the absolute value of spectra and x^ denotes the absolute value of post processed 
NMR spectra. For evaluation, T is set to two values. First, T is set to zero, which means a spectrum 
with small absolute values, possibly noise, are not suppressed. Second, T is set to the lowest value of 
contour when plotting the 2D spectrum. This is reasonable because peaks with absolute values smaller 
than T are not seen in the contour plot. 

Suppose x denotes the reconstructed spectrum from undersampled FID, relative li norm error 
(RLNE) is defined to measure the reconstruction error as: 



where x is the reconstructed spectrum from fully sampled FID and 0 < x, x T < 1. RLNE evaluates the 
normalized error presented in the reconstructed spectrum from undersampled FID. The lower the 
RLNE is, the better the reconstructed spectrum is consistent to the fully sampled spectrum. 

4.1. Reconstruction of the spectra 

The improvement by using the proposed method is verified from the less crowed l H- l H COSY 

1 13 

spectrum and more crowded H- C COSY spectrum. The sampling patterns of the two spectra are 
shown in Figure 8. 

Figure 9(c-h) show the reconstructed l H- l H COSY spectra corresponding to the sampling pattern 
in Figure 9(a) with a sampling rate of 0.20. With the l\ norm minimization, all the peaks are recovered 
successfully by using identity matrix (Figure 9(d,f)), while some peaks are lost by using wavelets 
(Figure 9(c,e)). 




(23) 



RLNE = 



(24) 
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Figure 8. Sampling pattern used in simulation, (a) Cartesian sampling pattern with 
sampling rate 0.20 for the 2D l H- l H COSY spectrum (N\ = 256 points) in Figure 4(a); and 
(b) Cartesian sampling pattern with sampling rate 0.25 for the 2D ^-^C COSY spectrum 
(N\ = 128 points) in Figure 5(a). 





(a) 



(b) 



Figure 9. CS reconstruction of a 2D H- H COSY spectrum using wavelet and identity 
matrix. (a,b) reconstructed spectra using fully sampled FID and undersampled FID with 
zero filling, respectively; (c,d) reconstructed spectra using wavelets and identity matrix 
with IST-based l\ norm, respectively; (e,f) reconstructed spectra using wavelets and 
identity matrix with PSOCA-based l\ norm, respectively; (g,h) reconstructed spectra using 
wavelets and identity matrix with PSOCA-based t p norm, respectively. 



^0.5 
f 1.5 
^2.5 
2 3.5 
« 4.5 
E 5.5 
£6.5 

ffi7.5 
- 8.5 

8 



^0.5 

f 1.5 
p. 

3 2.5 
I 3.5 

in 

15 4.5 



n r 



i i 



o ♦ 



•v. 



/ • 



i i i i i 



1 

0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 



5 7.5 6.5 5.5 4.5 3.5 2.51.5 0., 
*H Chemical shift (ppm) 

(a) 



■g 5.5 
£6.5 
K7.5h 
" 8.5 
8 



i i i i i i 



5 7.5 6.5 5.5 4.5 3.5 2.51.50. 
'H Chemical shift (ppm) 

(c) 



0.5 
f 1.5 
3 2.5 
i 3.5 

GO 

73 4.5 
£ 5.5 h 
£6.5 
^7.5 
" 8.5 



"i r— i 1 a i ! 1.41. 



4 j < 



V 

1 Ufa 



.5 7.5 6.5 5.5 4.5 3.5 2.51.5 
*H Chemical shift (ppm) 

(b) 



0.5 




.5 7.5 6.5 5.5 4.5 3.5 2.51.51 
'H Chemical shift (ppm) 

(d) 



Sensors 2011, 11 



8902 



Figure 9. Cont. 
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Since the contours for the marked peaks look faint, we also plot the ID slices along the indirect 
dimension in Figure 10. The height of one peak in the wavelet-based reconstruction in 
Figure 10(a,b) are much lower than those in the fully sampled spectrum, leading to the peak lost in the 
contour plots in Figure 9(c,e). 

Figure 10. ID slices along the indirect dimension for the chemical shift of 8.2 ppm (a-c) 
or 7.2 ppm (d) in the direct dimension, (a) Spectra reconstructed with IST-based l\ norm; 
(b) spectra reconstructed with PSOCA-based l\ norm; (c) spectra reconstructed with 
PSOCA-based ^0.5 norm; (d) spectra reconstructed with PSOCA-based ^0.5 norm. 
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Figure 10. Cont. 
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Furthermore, the nonlinear operation on wavelet coefficients induces the artifacts labeled in 
Figure 9(c,e). This phenomenon is also observed in the ID slices shown in Figure 10(a,b), where 
wavelet reconstruction generates illusive peaks. With the 4.5 norm minimization, the errors caused 
from wavelet and identity matrix reconstruction are reduced, as shown in Table 1. One can still 
observe the reduced peak height and artifacts in wavelet-based reconstruction, but identity matrix 
performs very well (Figure 10(d)). The advantage of 4.5 norm over l\ norm is obvious in the crowded 
^-^C COSY spectra, as will be shown in the following discussion. 

Table 1. Reconstruction error of a 'H-'H COSY spectrum. 



Methods 



Zero-filling 



1ST I 



PSOCA I 



PSOCA A).s 



Wavelet 



Identity matrix 



RLNE 
(T = 0) 
RLNE 

(T = 0.1) 
RLNE 
(T = 0) 
RLNE 

(T = 0.1) 



2.054 
0.059 
2.054 
0.059 



0.415 
0.012 
0.282 
0.010 



0.393 
0.010 
0.273 
0.007 



0.430 
0.007 
0.245 
0.022 



Figure 11 shows the reconstructed ! H- 13 C COSY spectra corresponding to the sampling pattern in 
Figure 8(b) with a sampling rate of 0.25. Some peaks are obviously lost in the reconstructed spectra 
using wavelets with both t\ norm and £0.5 norm minimization (Figure ll(c,e,g)). These lost peaks are 
found in the identity matrix-based reconstruction spectra (Figure ll(d,f,h)). With the £0.5 norm 
minimization, the intensities of the peaks marked with arrow in Figure 1 1(h) are more consistent to the 
fully sampled spectra in Figure 11(b) than those in the reconstructed spectra with the £\ norm 
minimization (Figure 1 l(d,f)). The smallest reconstruction error is achieved with the proposed identity 
matrix-based £0.5 norm minimization method (Table 2). 
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1 13 

Figure 11. CS reconstruction of a 2D H- C COSY spectrum using wavelet and identity 
matrix. (a,b) spectra reconstructed using fully sampled FID (N\ = 128 points) and 
undersampled FID with zero filling, respectively; (c,d) spectra reconstructed using 
wavelets and identity matrix with IST-based t\ norm, respectively; (e,f) spectra 
reconstructed using wavelets and identity matrix with PSOCA-based l\ norm, respectively; 
(g,h) spectra reconstructed using wavelets and identity matrix with PSOCA-based £0.5 
norm, respectively. 
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Figure 11. Cont. 
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Table 2. Reconstruction error of a 'H-^C COSY spectrum. 



Methods 



Zero-filling 1ST h PSOCA A PSOCA l 0 . 5 



Wavelet 



Identity matrix 



RLNE 
(T = 0) 
RLNE 

(T = 0.1) 
RLNE 
(T = 0) 
RLNE 

(T = 0.1) 



1.687 
0.098 
1.687 
0.098 



0.547 
0.044 
0.422 
0.033 



0.533 
0.042 
0.405 
0.031 



0.541 
0.042 
0.343 
0.027 



All above simulation results demonstrate that wavelet-based reconstruction obviously induces the 
loss of some peaks in the crowded 'H- 13 C COSY spectrum and loss of some weak peaks in the less 
crowded l H- l H COSY spectrum. The wavelet may even worsen the reconstructed spectra. Thus, it is 
not a good choice to use wavelets for the self-sparse spectra discussed in this paper. 

4.2. Discussion on the Computation 

Our simulation is run on a dual core 2.2 GHz CPU laptop with 3 GB RAM. The computational time 
for the algorithms using wavelet is two times that using the identity matrix, as shown in Table 3. 



Table 3. Running time for reconstruction of a NMR spectrum (unit: second). 



Methods 


Zero-filling 


1ST C\ 


PSOCA it 


PSOCA e 0 s 


'H^H ! h- 13 c 


^-'h 'h- 13 c 


l H- l H ! H- 13 C 


l H- l H ! H- 13 C 


Wavelet 


0.1 0.1 


11.1 56.8 


8.5 70.4 


29.1 221.2 


Identity matrix 


0.1 0.1 


5.9 27.5 


5.7 31.8 


16.0 105.6 



In the simulation, with the gradual increase of continuation parameter /?, the previous solution was 
used as a 'warm start' for the next alternating optimization in the PSOCA. For a given /?, with the 
increase of iterations in inner loop, the difference between reconstructed spectra decreases (see 
Figure 12(a)), so does the error between the reconstructed spectrum and the fully sampled spectrum 
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(see Figure 12(b)). The reconstruction error decreases when /? becomes large in the outer loop. The 
computational time of ^0.5 norm minimization in PSOCA is nearly four times as that of l\ norm 
minimization, as shown in Table 3. 

Figure 12. Numerical performance of PSOCA. (a) The £2 norm of difference between 
reconstructed spectra in the current and previous iteration when J3 = 2 12 in inner loop; 
(b) the reconstruction error RLNE of the reconstructed spectra when J3 = 2 12 in inner loop; 
and (c) the reconstruction error RLNE versus the iterations in outer loop in PSOCA. 
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5. Conclusions and Future Work 



Random sampling in the indirect dimension is introduced to reconstruct 2D self-sparse NMR spectra 
within the CS framework. Based on the assumption of sparsity of NMR spectra, one may remove the 
aliasing by penalizing the l\ norm on the coefficients of the sparse representation of NMR spectra. 
Considering the sparsity and the coherence property, we demonstrate that wavelet transform may 
reduce the peak height and result in loss of peaks. Thus, a wavelet is not necessary and even worsens 
the reconstruction of self-sparse NMR spectra. With the £ p (p = 0.5) norm minimization, the quality of 
reconstructed spectra can be further improved. 

However, how to define the meaningless peaks depends on applications and a qualitative analysis of 
self-sparse NMR spectra is needed in order to satisfy the requirement of CS. By defining regularity of 
ideal Lorentizian peaks with aspect to typical vanishing moment wavelet basis, it is possible to give a 



Sensors 2011, 11 



8907 



boundary for the approximation error of Lorentizian peaks in wavelet representation. Thus, one may 
quantify the sparsity of spectra composed of ideal Lorentizian peaks using wavelets. Another way is to 
set up a database and analyze the sparsity of the meaningful peaks based on the prior knowledge of 
chemists. Since the peak height may be reduced in the wavelet-based reconstruction and this reduction 
depends on the crowd of peaks, it is expected to give a quantitative analysis on the effect of 
using/skipping wavelet transform by setting up a simulated spectrum or spectrum from real chemical 
substance, in which the crowd of peaks and the fixed relative height of peaks are pre-defined in 
the spectrum. Besides, based on the coherence property in CS, the analysis of the performance of 
different random sampling schemes, e.g., Poisson disk sampling, may lead to further reduction of 
sampling rate and reconstruction error. Extension of the proposed method on higher dimensional NMR 
spectra is worth investigating. 
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